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Abstract 

The statistical properties of the local optima (metastable states) of the infinite 
range Ising spin glass with p-spin interactions in the presence of an external mag- 
netic field h are investigated analytically. The average number of optima as well 
as the typical overlap between pairs of identical optima are calculated for general 
p. Similarly to the thermodynamic order parameter, for p > 2 and small h the 
typical overlap qt is a discontinuous function of the energy. The size of the jump in 
qt increases with p and decreases with h, vanishing at finite values of the magnetic 
field. 



Short Title: landscape of the p-spin model 
PACS: 05.50+q, 87.10-He, 64.60Cn 
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1 Introduction 



The emphasis professed by KaufFman on the role of the topology of the fitness 
landscape as a source of order in contraposition to natural selection has arisen 
considerable interest in the study of the statistical properties of fitness landscapes 
[0. The central issue is the limitation imposed by the structure of the fitness 
landscapes on adaptive evolution, viewed as a local hill-climbing procedure via fitter 
mutants. (See for a lucid criticism of these ideas.) For sake of concreteness, let us 
consider a population of asexually reproducing haploid organisms whose genotypes 
are described by sequences of N Ising spins s = (si, . . . , s^v) with Si = ±1. In the 
discrete space of the 2^ possible sequences, evolution is modelled by an adaptive 
walk defined as a connected walk through a succession of neighboring sequences 
(i.e., sequences that differ by a single spin only) each of which possessing improved 
fitness . There are several questions of interest whose answers may shed light on 
the structure of the landscapes as, for instance, the number of fitness optima in the 
sequence space and the similarity between these optima, to mention only those we 
will address in this paper. 

Most of the analyses have concentrated on the NK model of random epistatic 
interactions since it possesses a tunable control parameter that regulates the rugged- 
ness of the fitness landscape [0, An alternative (and more appealing to the 

physicists) class of fitness functions was proposed by Amitrano et al. ^ , namely, 
the Ising spin glass with p-spin interactions defined by the random energy function 
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where the coupling strengths are statistically independent random variables with a 
Gaussian distribution 
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and h is the external magnetic field. In this context the fitness value ascribed to 
a sequence or genotype s is the reverse of its energy. Thus the fitness maxima 
correspond to the energy minima of (|l|). Henceforth we will refer to the fitness 
maxima or energy minima as simply optima. For p ~ 1 01 h —^ 00 the energy 
(P gives a single-peaked, smooth correlated landscape, while the limit p ^ 00 
corresponds to the random energy model of Derrida and yields an extremely 
rugged, uncorrelated landscape. The case p = 2 is the well-known SK model 
which exhibits a large number of highly correlated local optima |l^ . 

For general p, little is known about the statistical features of the landscape 
generated by the energy function (|^). A result worth mentioning is that, for h = 0, 
the correlation between values of Hp for different configurations is given by |l^ 



(3) 
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where 

1 ^ 

is the overlap between the two arbitrary states s° and s**. Here the average indicated 
by (. . .) is taken over the probabihty distribution of the couphngs (^). Thus, as 
mentioned before, the correlations between energy levels vanish for p — > oo. 

The thermodynamics of the p-spin Ising model has been investigated within the 
replica framework |l^, In particular, for p = 2 the order parameter function 
q{x) tends to zero continuously as the temperature approaches a critical value at 
which the transition between the spin glass and the high temperature (disordered) 
phases takes place 0. For p ^ oo, the system has a critical temperature Tc 
at which it freezes completely into the ground state: q{x) is a step function with 
values zero and one, and with a break point at x — T/Tc 0. The situation for finite 
p > 2 is considerably more complicated. There is a transition from the disordered 
phase to a partially frozen phase characterized by a step function q(x) with values 
zero and qi < 1. As the temperature is lowered further, a second transition occurs, 
leading to a phase described by a continuous order parameter function ]T^ , [Tst . 

The goal of this paper is to investigate the statistical properties of the fixed 
points (local or global optima) of adaptive walks on the fitness landscape defined 
by equation (|^). The energy cost of flipping the spin Si is STip = 2Ai where 

^ ^ Jii^.-.ip ^i'^i2 • • ■ ^'ip ~^ hSi- (5) 
i^2<■■■<ip 

is termed the stability of Si. Since in an adaptive walk only flippings or moves that 
decrease the energy (i.e., increase the fitness) are allowed, any state s that satisfies 

A, > Vi (6) 

is an optima of the fitness landscape. Clearly, counting the number of states that 
obey (^ is equivalent to calculating the number of solutions of the zero-temperature 
limit of the celebrated TAP equations iQ. For non-zero temperature, the quite 
involved calculation of the average number of solutions of the TAP equations has 
been carried out for p = 2 [|l^ as well as for general p |]l6|. However, systematic 
analyses of the typical energy of the local optima and of the effects of the external 
magnetic field have been undertaken for the simplest case only, namely, p = 2 at 



zero temperature |17, 1^, We note that in the statistical mechanics context the 
local optima are usually termed metastable states. 

In this paper we study at length the effects of the magnetic field h on the 
structure of the local optima of the p-spin energy landscape. More pointedly, we 
calculate analytically the average number of local optima with a fixed energy density 
e, denoted by {Af{e)). Although this analysis is quite straighforward, it is justified 
since the dependence of that quantity on e and h has not been investigated for 
general p. In fact, we note that results of extensive numerical simulations aimed at 
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measuring (A/'(e) have been reported recently |Q . More importantly, we calculate 
the average number of pairs of local optima with overlap q and fixed energy density 
e. This quantity, denoted by {M{q,e)), allows us to determine the typical overlap 
Qt between pairs of local optima with energy density e. Since {M{qt,e)) is directly 
related to the second moment of Af{e) , we can determine the regions in the space 
of parameters (p, e, h) where this random variable is self-averaging. 

The remainder of the paper is organized as follows. In Sec. we derive the 
formal equation for the nth moment of the random variable A/'(e). Then we use 
that result to calculate the average number of local optima (A/'(e)) in Sec. ||, and 
the average number of pairs of local optima {A4{q,e)) in Sec. |j. Finally, some 
concluding remarks are presented in Sec. ^. 



2 The formalism 

The number of local optima M{e) with fixed energy density e can be calculated by 
introducing the quantity Yg defined by 

Y - f 1 iieN = np (s) and A, > Vi , 
^ 1^ otherwise ' 

so that 

AA(e) = Tr^Fs, (8) 

where Trg denotes the summation over the 2^ states of the system. We are inter- 
ested in the evaluation of the moment ([^(e)]") for n — 1,2, which can be written 
as 

= Tr,i...Trs.W(y,i =l,...,ys" =1) (9) 

where W (Ygi = 1, . . . , Ys'^ = 1) is the joint probability that the n random variables 
Ygi , . . . , Ygi assume the value 1. Using the definition 

w = 1, . . . , y.. = 1) = / n <5 [e - (s'^) /N] II e (A^ ) , (lO) 

\a=l i I 

the equation for the nth moment becomes 

{W{^)T) - ( n Tr,.<5 [eN - Hp {s^)] [] 6 (A,-)\ . (11) 

\a=l i I 

where Q{x) = 1 if a; > and otherwise. We have presented the derivation of 
equation (fill) in detail because some authors have written the random variable 
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A/'(e) in terms of the delta function directly jl^, Clearly, this procedure is 

correct only for the moments of Af{e) as shown above. 

In the next two sections we concentrate on the explicit evaluation of equation 
( pi] ) for n = 1 and 2. To facilitate those calculations, we express the energy Tip (s) 
in terms of the stabilities A^, 

nJ,{s)^--Y,i^^ + hip-l)s^), (12) 

^ i 

SO that the dependence on the couplings in equation (^) appears only through the 
stabilities A^. 



3 Average number of optima 



Using the integral representation of the delta function and the auxiliary relation 
(^2!) we can write the first moment of A/'(e) as 



m^)) = |Jexp(iiV«)ny'°° ^^e(A,)exp(iA,A.) 

xTr^exp -ih^AiS, + -e^ (A^ + h{p - l)si) 

i ^ i 

x(exp -i^A; ^ J^^2...^pS^S^2...S,^ 
\ y i i2<---<ip 

The average over the couplings can be easily carried out using the identity 

/ P \ 



^ ^ ^ ^ Jii2...ipSiSi2 ■ ■ 

i i2<---<ip 

and yields, in the limit N ^ 00, 



ii<...<ip \fe=l / 



exp 



exp 



ii<...<ip \fc=l / 



(13) 



(14) 



(15) 



The remaining calculations are straightforward: a Gaussian transformation allows 
us to decouple the sites in (|l5|), so that the integrals over A, and A, as well as the 
trace over the spins can be readily performed. As usual, we conclude the calculation 
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by carrying out a saddle-point integration over two appropriately rescaled saddle- 
point parameters. The final result for the exponent / in (A/'(e)) = e^^ is 



In 



c'^'^erfc (-/X - /i) + c-'"'erfc (-/x + h) , (16) 



where 



(17) 

Here the saddle-point parameters v and /i are obtained by solving the equations 
df /dv = and df/dii — simultaneously. In figure 1 we present the exponent / 
as a function of e for p ^ 2 and several values of h. For sake of clarity we present 
only positive values of /. The decrease in the number of local optima as h increases 
indicates that the landscape becomes smoother, as expected. The results for p > 2 
are qualitatively similar, except that the peaks are higher and slightly broader. Two 
values of the energy density are particularly important, namely, the value at which 
/ reaches its maximum value ft, denoted by ef, and the lowest value of e for which 
/ vanishes, denoted by eq. While et gives the typical value of the energy density 
of the local optima, eo gives a lower bound to the ground state energy density of 
the spin model defined by the hamiltonian (Q) jl^. In figures 2 an 3 we present 
et and ft , respectively, as a function of h for several values of p. These quantities 



are easily obtained by setting = in equation (16). The single saddle-point 
equation df /dfi = possesses either one root (for either small or large values of 
h) or three roots (for intermediate values of h). The discontinuity in et that can 
be observed in figure 2 for p > 7 is due to the simultaneous disappearance of two 
of those roots. For p — > oo and finite h we find et — > (Hp) = 0, signaling thus the 
emergence of the so-called complexity catastrophe, i.e., the energy density of typical 
local optima equals the expected energy of a randomly chosen state [Q . We note 
that {Af{et)) — exp(/f A^) yields the average number of optima regardless of their 
energy values, i.e., the same result is obtained by dropping the energy constraint 
in the definition of Yg given in equation (Q) . In figure 4 we present eg as a function 
of h for several values of p. Clearly, since in the limit /i — > oo there is only one 
optimum, namely, s = 1, we find eo — > et = ^h. It is important to note that for 
p ^ oo, eo tends to a non-zero limiting value. This result illustrates the fact that 
the complexity catastrophe phenomenon affects the typical optima only. In fact, 
the increase of p has little effect on the ground-state lower bound eo, which for 
h = decreases from -0.791 for p = 2 ^ to - Vb[2 « -0.832 for p ^ oo §. 
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4 Average number of pairs of optima 

We define the numl 
energy density e as 



We define the number of pairs of optima with overlap q = —1,-1 + -^,...,1 and 



ATI ■ • • J ■ 



(18) 



where 6{m, n) is the Kronecker delta and is given by equation (|^). Following the 
procedure presented in Sec. ||, the average of M over the couplings is cast into the 
form 

{M{q, e)) = \ ^Tr,.Tr,.5 ^q, ^ s\s^ H b [eiV - (s'')] H 6 (A^)^ (19) 

The integral representations of the delta function and the Kronecker delta allow us 
to write this equation as 



(X(<Z,e)) 



2 y , 2! ^'^'^^^ n y ^ ^ exp [INe'^n 
X n Trs» ^^^e (An exp (iA^ A^) 

-\qY, s\sj -ihJ2 A^- + ^ 5] (A^ + % - l)sf) 



X exp 



X ( exp 



(20) 



ai i2<...<*p 



As in the previous section, the average can be performed with the aid of an identity 
analogous to (p^, yielding 



(...)= exp 



T E 

ii<...<ip 



E Ea^ 4 



.a=l \k=l 



(21) 



After some algebra, the argument of this exponential is rewritten in the limit N 

00 as 



{...} 



E 



E(a^ 



N 



E^^ 



pq 



p-i 



p{p-l)qP-^ 



2N 



(22) 
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The next step is to introduce via delta functions the auxihary parameters: Nmi = 
Ai, Nm2 = E» Af , Nvi = J2, ^}s}sl Nv2 = 'E^ Af s^s^, and their respective 
Lagrange multiphers in order to decouple the variables sf and Af for different 
sites i. Then the integrals over Af and A", and the trace over sf can be easily 
performed. As before, the auxiliary parameters as well as the Lagrange multipliers q 
and e are integrated out via a saddle-point integration. This part of the calculation 
is straightforward and quite unilluminating so we do not present any further detail. 
To proceed further we assume that the symmetry <-!■ between the two replicas 
remains intact, i.e., mi = m2 and vi = V2- This is a quite sensible assumption since 
the breaking of the replica symmetry that pervades the thermodynamic calculations 
[0, |l^, |l^ is very probably a consequence of the limit where the number of replicas 
goes to zero. In any event we will, conservatively, restrict the forthcoming analysis 
to pairs of identical optima only. The final result for the exponent g in {Ai{q, e)) = 
^ exp((/7V) is written more simply in terms of a new set of saddle-point parameters 
that are linear combinations of those introduced above. We find 



1 



VP 



qz 



2{p-l) 



(x + yy + q^-P (x^yY + il + qP)- 

4p 



2(p-i: 



[(1 + g) a; + (1 - (?) y] + In S {ly, x, y, z) - In 2 



(23) 



where 



s = e 



uh- 



Dt erfc 



— x — h 



x + h + qP-H 



+e" 



+e^ 



+e^ 



Dt erfc 



x+h 



Dt erfc 



-V+h 



Dteiic 



y-h 



.2p-2 



y + h-qP-H 
y/1 - q^P~^~ 



iP-H 



2p-2 



(24) 



Here Dt = dte * /-/tt is the Gaussian measure and h is given by (0). The saddle- 
point parameters v, x, y, z must be determined so as to maximize g. This is achieved 
by solving the four coupled saddle-point equations dgjdv — 0, dg/dx — 0, dg/dy = 
0, and dg/dz = 0. For g = 1 we find y = Q and hence g — f, a,s expected. 
Furthermore, for g = and h = we find x = y and z = so that g = 2/. Once 
{M{q, e)) is known, the second moment of A/'(e) can be calculated using the identity 



{M{qt,e)), 



(25) 
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since the sum is dominated by the overlap q — qt that maximizes equation ( p3| ) in 
the hmit N ^ oo. Hence we have {[J^{e)f) = exp (f'-'^^N) with given by (|3|) 
calculated at qt- Thus for h — the variance of the random variable Af{e) vanishes 
in the thermodynamic limit, provided that qt = 0. We note that although g = is 
always a point of maximum of g for h = 0, that maximum may not be the global 
one and, in that case, qt ^ 0. 

For fixed q, the dependence of 5 on e is similar to that shown in figure 1. 
Likewise, the maximum of g with respect to e, denoted by gt-, is determined by 
setting v = 0. In figure 5 we show this maximum as a function of q for p — 7 and 
several values of h. The quantity ^ exp (gt-ZV) can be viewed as the number of pairs 
of identical optima (in the sense that their energies and saddle-point parameters 
are identical) with overlap q, regardless of the specific value of their energies. For h 
not too large there appears a minimum for g « 1, indicating that around a typical 
optimum there is a region where other optima are rarer. The picture that emerges 
is one of clusters of many optima surrounded by comparatively smoother valleys. 
The typical energy ej of these optima is shown in figure 6 as a function of the 
overlap q. The typical overlap qt between the optima increases from zero at ft. = 
to one in the limit h ^ 00 since, as expected, the external magnetic field induces 
correlations between the optima. This is shown in figure 7, where we present qt as 
a function of h for several values of p. The discontinuity that appears for p > 7 is 
caused by the competition between the two maxima shown in figure 5. 

Next we consider the dependence of the typical overlap between identical optima 
on their energies. This analysis is more involved since, besides the four saddle-point 
equations, we have to solve the equation dg/dq = too. In figures 8 and 9 we show 
qt as a function of e for p — 2 and p = 3, respectively, and several values of the 
external magnetic field. For h = 0, in both cases we find = up to a certain 
value of the energy density (e = —0.672 for p = 2 and e = —0.792 for p = 3). Thus, 
as mentioned before, A/'(e) is self-averaging in this regime. Our results for p = 2 
are remarkably similar to those found in the replica calculation of the quenched 
average (In A/"(e)), with the typical overlap qt replaced by the saddle-point parameter 
q = ((si)g) |1^. Here (. . .)<: means an average over optima with energy density e. 
In particular, q vanishes for e > —0.672, indicating thus that A^(e) is self-averaging 
in this regime, in agreement with our results. However, while for p = 2, qt increases 
continuously from zero, for p = 3 there is a discontinuity at e = —0.792. The 
same phenomenon is observed for p > 3, with the size of the jump in qt increasing 
with p. This finding is reminiscent of the jump in the order parameter found in 
the thermodynamic calculations for p > 2 ^ . The discontinuity in qt can be 
understood by studying the dependence of the exponent g on the overlap q for p = 3 
and h = 0, shown in figure 10. Since the typical overlap is associated to the global 
maximum of g, the competition between the maximum at g = and the maximum 
at g > originates the jump in qt, which takes place at the energy density where the 
two maxima have precisely the same height. The situation for non-zero h is more 
complicated. The correlations induced by the magnetic field destroy the region 
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of self-averageness of A/'(e) . Interestingly, for a given h > there is value of the 
energy density for which the typical overlap is minimal. For p — 3 the effect of the 
magnetic field is to decrease the size of the jump in qt till it disappears altogether 
for h « 0.29. The results for p > 3 are qualitatively similar to those for p = 3. We 
mention only that the larger p, the larger the value of e at which the discontinuity 
occurs, and the larger the value of h at which it disappears. Unfortunately, the 
enormous difficulty of solving the system of five coupled equations prevents a more 
systematic analysis of these discontinuities. 



5 Conclusion 

The analytical investigation of the statistical structure of the energy landscape of 
the p-spin Ising model presented in this paper is of interest from the viewpoint of 
the traditional statistical mechanics of disordered systems |l^ as well as 

from the perspective of the study of adaptive walks in rugged fitness landscapes 
|ll|, 20 . Besides extending the calculation of the average number of optima to 



general p and non-zero magnetic field, we have focused on the characterization of the 
typical overlap qt between pairs of identical optima. Interestingly, the dependence 
of qt on the energy density e is reminiscent of the dependence of the thermodynamic 
order parameter on the temperature T ||l^, |l^ . We must note, however, that there is 
no relation between T and e since lnA/'(e) is not the entropy of the spin system. The 
quite complex effect of the magnetic field on the statistical properties of the energy 
optima motivates a more detailed study of the thermodynamics of the p-spin model 
for non-zero h. In fact, even the unambitious analysis of the first moment (A/'(e)) 
has unveiled an interesting interplay between h and p that lead to a discontinuity 
in the typical energy density of the optima. Moreover, we have found that the 
magnetic field decreases the size of the jump in the typical overlap qt that occurs 
for p > 2. It would be interesting to investigate whether a similar effect occurs 
for the thermodynamic order parameter as well, which might lead, eventually, to a 
continuous phase transition. 

To conclude, we must mention that the calculations presented in this paper are 
free of all the mathematical subtleties that permeate the replica analyses of the 
infinite range Ising spin glass [|[ |o) . Thus our results present a reliable account of 
the statistical properties of the p-spin energy landscape which, though may have 
little relevance to the thermodynamics of the model, are of considerable interest 
to the characterization of the fixed points (metastable states) of adaptive walks 
(zero-temperature Monte Carlo dynamics) on that landscape. 
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Figure captions 



Fig. 1 The exponent / in {Af{e)) = e-^^ as a function of the energy density e for 
p = 2 and /i = 0, 0.5, 1.0, and 1.5. 

Fig. 2 The typical energy density et of the local optima as a function of h for (from 
top to bottom) p = 2 to p = 10. For p — > oo we find et — > 0. The dashed straight 
line is Ct = —h. 

Fig. 3 The exponent ft in the expression for the average number of optima 
{Af{et)) = c-'*^ as a function of h for (from bottom to top) p = 2 to p = 10. 
For p oo we find ft In 2. 

Fig. 4 The lower bound eo to the ground state energy density as a function of h 
for (from bottom to top) p = 2, 3, 4, and oo. The dashed straight line is eo = —h. 

Fig. 5 The exponent gt in the expression for the average number of pairs of 
identical optima {M{€t, q)) = e^*^ as a function of for p = 7 and (from top to 
bottom) h = 0,l, 2, 2.5, 3, 3.3, 3.6, 3.8, 4 and 4.2. 

Fig. 6 The typical value of the energy density of a pair of identical optima as a 
function of the overlap g for p = 7 and (from bottom to top) h = 0, 1, 2, 2.5, 3, 
3.3, 3.6, 3.8, 4 and 4.2. 

Fig. 7 The typical value of the overlap between pair of identical optima as a 
function of h for (from left to right) p = 2 to p = 8. 

Fig. 8 The typical value of the overlap between pair of identical optima as a 
function of their energy density for p = 2 and h = 0, 0.5, 1.0, and 1.5. The marked 
points correspond to /^^^ = 0. 

Fig. 9 Same as figure 8 but for p = 3, and ft = 0, 0.27, 0.5, 1.0, and 1.5. 

Fig. 10 The exponent g in the expression for the average number of pairs of 
identical optima {A4{e^q)) = e^^ as a function of g for p = 3, ft = 0, and (from top 
to bottom) e = -0.73, -0.75, -0.77, -0.79, and -0.81. 
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